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0. Introduction 


During the last decade, discrete methods other than the classical 
finite-differences — have gained an increasing popularity while used for the 
approximate solution of time-dependent problems* Most noticable, are the 
(pseudo) spectral and Galerkln methods, e.g. [4-9], [13-14], [19], [21] and the 
references therein • 

The purpose of these notes is to give a unified survey on these three 
classes of discrete methods — finite differences, spectral and Galerkin, 
discussing some of the theoretical aspects with regard to their accuracy, 
stability and efficient implementation* As a model problem for our 
discussion, we consider the one-dimensional S3rmmetrlc hyperbolic system 

(0.1) 3^u(x,t) = A(x,t) 3^u(x,t) + B(x,t)u(x,t) , A(x,t) - A^(x,t)* 

Here and elswhere in the paper, w^ denotes the transpose of a given vector, 

* . * I/9 

w its conjugate transpose, and Hwll=(w w) ^ its Euclidean norm; similar 

notations are used for matrices. We will also briefly mention systems of 

parabolic type, which are more favored by the kind of arguments discussed 

below, due to the presence of dissipation* To avoid further complications 

that arise with time-discretization and handling boundary conditions, we 

restrict our attention to the periodic method of lines* This allows us to 

make a rather detailed study of various types of approximations to (0.1), 

where the spatial dif ferentation is replaced by its discrete counterpart. 

We begin, in Part I, discussing finite difference methods* Our approach 
— slightly different than usual — emphasizes the matrix representation of 
such methods* The reason is a two-fold one: first, the standard approach via 

Von Neumann analysis is by now classical and can be found in a variety of 
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references , e*g. [15] ; second, viewing these discretizations in the language 
of matrices allow us to move quite naturally to our discussion on spectral and 
Galerkin methods in the second and third parts. Indeed, the generality of the 
abstract discrete dif ferentation operator dealt with in Part I will prevail 
for finite-difference as well as spectral and Galerkin methods. One of the 
main ob j ec tives of this review is , in f ac t , to show the intimate relation 
between the three: the spectral-Fourier method can be viewed as a special 
centered finite differencing based on an ever increasing number of gridpoints 
peridoically extended, and both result from an appropriate choice of basis 
functions used in the Galerkin method. 

We focus our attention on the all important question of stability. It is 
shown that antisymmetry periodicity as well as a locality restriction are 
essential properties that a discrete differencing method should share with the 
differential problem, for the resulting discrete system to be stable. The 
locality restriction can be equivalently expressed by the boundedness of 
amplification blocks associated with the highest modes. The accuracy 
requirement, on the other hand, is determined by the exactness of differencing 
the lower modes. The combination of the two guarantee convergence, as the 
lower modes carrying most of the information are accurately represented, while 
the highest modes are not, yet stability assures us that they are not 
amplified and hence rapidly tend to zero, just as is the case with the 
differential problem. 

Both properties of accuracy and stability are well accommodated in 
discrete methods having finite degree of accuracy; they contradict each other, 
however, with highly accurate methods. We are then led to a discussion on the 
skew-symmetric differencing and smoothing procedures. Both aim at dissolving 
this contradiction by bounding the highest modes' amplification blocks , yet 
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leaving the lower, accurate modes, unharmed. 

In Part II we amplify these points with regard to the spectral-Fourier 
method, from still a slightly different point of view. We base the whole 
stability analysis in this case on the aliasing formula, relating the Fourier 
coefficients of a given periodic function to those of its equidistant 
interpolant. This is the single most important formula, which dominates the 
question of accuracy versus stability in this case. It naturally arises with 
the Fourier method as the aliasing dilemma, and its usual remedy is again by 
either skew-symmetric differencing or via smoothing. From this point of view, 
finite difference methods having finite degree of accuracy can be viewed as 
special cases of the Fourier method with a built-in smoothing which guarantees 
their stability. 

In Part III we discuss Galerkln-type methods. Again there is an emphasis 
on the close connection with finite-difference and spectral-Fourier methods. 
Stability follows in this case, due to lack of aliasing . Once the exact 
Fourier coefficients are discretized, we find ourselves dealing with exactly 
the same kind of arguments introduced before. 

Finally, in order to make these notes self-contained, we have collected 
in the Appendix some basic properties of Toeplitz and circulant matrices; 
these play a vital role in the foregoing analysis. 
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Part I. Finite Difference Methods 


1. Finite Difference Operators 

Let v(x) be a 2ir-perlodlc m-dlraenslonal vector function, whose values 
= v(x^) are assumed known at the gridpoints = vh, h = — 

V = A second order accurate approximation to its 

derivative, D^v(x), is given by the centered divided difference 

When augmented by the periodicity of v, these divided differences are well 
defined at all gridpoints x = x^, v = 0,1,***,N-1. The transformation which 
takes the vector of the assumed known gridvalues v = [vq,***,v^ into the 

vector of divided differences [y] = [d^ ( h) [ vq] , • • • ( h ) )' is 

linear, and hence has a matrix representation 
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Its entries are being la-dimensional blocks. Similarly, fourth and sixth order 
accurate centered divided differences are given by 


4D (h) - D (2h) 

D^(h)[v(x)] ^ tv(x)] 


(I.I4) 


8[y(xfh)-v(x-h)] - fv(x+2h)-v(x-2h) 1 
12h 





accurate centered divided difference given by [9, Section 3] 


(l*l2g) 


" (s+k) !(s-kjl’ 


k-1 


likewise. It has an antisymmetric block clrculant matrix representation 
^28 " 


( 1 . 22 ^) 


9_ [v] =* D* V. 

zs 


As s Increases, so does the amount of work required to perform the 
multiplication on the right-hand side of (1. Traditionally , f Inlt e- 
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a total amount of work of N*s operations (1 operation = vector addition + 
vector multiplication by a scalar). For large values of s, becomes a 

full matrix whose multiplication requires an increasing amount of work, up 
to N operations. This number of operations can be reduced, however, by 
taking into acount that the matrix is a circulant one, and as such, can 

be diagonalized by the block Fourier matrix F. Specifically, let the block 
Fourier matrix F be given by 

(1.4a) [f]j^ = 0 < j,k < N-1 

with the conventional notation 

(1.4b) = ^-n, n = integral part of N/2, 


to be used throughout the paper. We then have (see the appendix for details) 


( 1 • 3 ^ 22 ) 


D. = NF A. F 
'^2s ^2s 


with the block diagonal matrix = ^2s^^^ given by 


a-3b2^ (ijsljj 





Multiplication of in its spectral representation (^*32g) can be 

efficiently implemented by two FFT's and N scalar multiplications which 
amount to SNlogN operations. 

In general, we consider an abstract discrete dif ferentation operator, 
whose matrix representation D = 2(h) is only required for the two basic 
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propertles of being an antisymmetric block circulant one 

(1.5) [sijk - 

(The antisymmetry requirement which corresponds to centered type differencing 
is in fact not as essential, but will suffice for our discussion below.) Such 
matrices admit the spectral representation 

(1.6) D = NF*AF 


with a block diagonal matrix A, whose diagonal consists of the so called 
amplification blocks 


(1.7a) 


[A] 


N-1 


ii 


E I cL , 

m , n K m 

k=0 


0<j<N-l; 


since D is assumed to be antisymmetric, d^ + “ 0, and hence 


(1.7b) 




n 


2i. I <L sin(j^kh), 
k=l 


0<j<N-l. 


The discrete dlf ferentatlon as given by the spectral representation of 
D, D“NF*AF may be Interpreted now as follows: from the grldvalues 

^v|0<v<N 1* have the discrete Fourier modes given by 


( 1 - 8 )* 




1 -iw^vh 

N ^ ® V, 

v=0 




then, each one of these modes is discretely differentiated as it multiplied by 
the ampllciation factor \ and finally, the differentiated modes 
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^ transformed back into the physical gridspace upon 

* _ -1 

multiplcation by NF = F 


2. Stability of Finite Difference Approximations 

Replacing the spatial derivative in (0.1) by its discrete counterpart 

D = D(h), we end up with the finite-difference approximation^^^ 

(2.1a) 3^v^(t) = A(Xy)D(h) [v^(t) ] + B(Xy)Vy(t); 

introducing the block diagonal matrices A = diag[ A(Xq) ,* • * ,A(x^ j^)] 

B = diag[B(xQ),***,B(x^_j^)], it can be put in a matrix form 

(2.1b) 9^v(t) = A3^(t) + By(t). 

The time-dependent difference equation (2.1) serves as an approximation to the 
differential problem (0. 1) , in the sense that any smooth solution, u, of 
(0.1), satisfies (2.1) modulo a small local truncation error ^(h) = x(h;t) 

(2.2) 3^u(t) = ADu(t) 4 * Bu(t) + T(h;t). 

The approximation is said to be accurate of order a if !lT(h)H =^[h^]. 
With D = for example, one obtains a difference approximation which is 

accurate of order 2s, To link the local order of accuracy 


^Also termed as the method of lines , to distinguish from the fully 
discretized problem in both space and time* 
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wlth the desired global convergence rate of the approximation, one has to 
verify Its stability. That Is, the approximation (2.1) Is said to be stable 
If for all sufficiently small h we have 

(2.3) Hexp[^t]ll < K = K^, 0<t<T. 

Observe that the stability definition is independent of the lower order 
terra. By, the reason being that stability in the above sense is, in fact, 
insensitive to such low-order perturbations. This is the content of the 
following classical perturbation lemma, whose proof is given here for 
completeness as it will play an essential role in our discussion below, (see 
e.g. [15, Section 3.9], [16] and under a much more general setup [20]). 

Perturbation Lemma Let A be a given linear operator such that 

llexp[At]ll < K^, 0<t<T. 

Then, altering A by adding a "low-order” bounded perturbation B, retains 
the exponent boundedness 


IIexp[(A4B)t]ll < K(t), K(t) = K^e ^ , 0<t<T. 


Proof The solution of the Inhomogeneous linear differential equation 


(2.4a) 


3^w(t) = Lw(t) +G(t), w(t=0)=w(0) 


Is given by 
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(2.Ab) w(t) = e^*^w(0) + / 

5=0 

Applying (2.4b) with L = A + B and G = 0, we then get 

w(t) = 

hence (2.4a) can be also written in this case as the inhomogeneous problem 
3^w(t) = Aw(t) + Be w(0) . Applying (2.4b) once more, this time with 
L = A and G = Be ^ w(0 ) , we obtain 

w(t) = el-^*^^w(0) + / 

5=0 

Equating the last two representations of w(t) which are valid for arbitrary 
initial data w(0), we arrive at the well-known identity 

t 

exp[(A+B)t] = exp [At] + / exp [A(t-5) ] *B*exp [ (A+B)5 ]d5 • 

?=0 

Taking norms on both sides we find 


t 

K(t) < K™ + K ‘IIBII* / K(5)d5; 

5=0 


by Gronwall 
and hence 


Inequality, we conclude that / K(5)d5 < UBN • [e ~l]» 

5=0 

, ^ K_* IIBII -t ^ K,_‘«B0*t 

K(t) < + K^- IIBII- IIBII" *[e -1 ] = K^e ^ 


which completes the proof. (We remark that similar arguments apply for the 
analogous question vdiich concerns the discrete framework discussed in [15, 
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Section 3.9]). 

From the preturbation lemma we see that stability is equivalent to the 
boundedness of or — what amounts to the same thing — to the 

continuous dependence of the discrete solution v(t) on its initial data 
X(0) 

[(Ag+B)t] 

(2.5) Ily(t)H = lie < K(t) • Hy(0) II ; 


indeed if stability holds in the sense that Hexp[^t]ll is bounded, see 
(2.3), then by the perturbation lemma with A = AD and B = B, so is 
Ilexp[(AD^+B)t] II. On the other hand, if Hexp [ (4E+S) ^ ] H is bounded then by the 
perturbation lemma with A = ^+B and B = -B, so is Hoxpf^St ] II . 

Granted stability, we can now estimate the global error 
E(t) = u(t) - y(t): subtracting (2.1) from (2.2), we find that it is governed 

by 


\E(t) = (^H-B)E(t) + T(h;t) 


and hence is of the form 


[(j^+B)t] t [(AD+B)(t-0] 

E(t) = e E(t=0) + / e T(h;5)d5; 

C=0 


using the perturbation lemma, we end up with the error estimate 

t 

llE(t)ll < K(t)*ilE(t=0)ll + sup llT(h;5)ll* / K(5)d?. 

0<5<t 5=0 

Thus, if an ot-order accurate approximation is initialized with ot-order 
accurate data, HE(t=0)H =<^[h^], its stability will retain ot-order of 
convergence later on, llE(t)H =<^[h^]. 
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Verifylng stability Is our main objective in the rest of this section* 
We start by rewriting 


ADt = V2 ()^B+RA)t +I/2 

where by the symmetry of A and antisymmetry of D, the first term on the 
right is antisymmetric and, therefore, has a bounded exponent 

llexp[ I /2 (^+DA)tlll = 1; 

hence, by the perturbation lemma, with the second term on the right viewed as 
a low-order perturbation of the first, the exponent of the sum of the two 
terms is bounded provided the second is 

II V 2 (^-DA) “ *T 

(2.6) llexp[^t]ll < K^e , 0<t<T. 

Thus we are left with finding a bound for the symmetric part of 

whose (p,q) block entry is given by 

[ I /2 (i^-DA)]pq =^/ 2 dfp_q] *[A(Xp)-A(x^)l, 0<p,q<N-l; 

since A(x) Is assumed symmetric 27T-periodic, H A(x^ )-A(x^) H 

< h* Max IIA^ (x) 11 •Mln[ Ip-q I ,N- Ip-q 1 ] , and hence V 2 (^*5^) is bounded 

0<x<2tt 

entrywise and therefore in norm, by the matrix whose (p,q) block entry is 
given by 

Y’MaxllA' (x) II • I *Mln[ |p-q| ,N-Ip-q 1 ] 
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This last matrix is a clrculant one. In the appendix we show, see Corollary 

(A. 8), that the norm of such matrix does not exceed (and In fact equals. In 

our case) the absolute value sum of Its elements along Its first (p = 0) 

row, — •MaxIlA' (x) II • 1 Mln[q,N-q ] • ] d |; recalling that D Is antlS 3 mimetrlc , 
q=0 ^ 

%-k “ finally end up with the desired hound 


n 

(2.7) "V2(^-5A)II < h* I kid I* Max llA'(x)H. 

k=l ^ 0 <x<2tt 


Insterting the last bound into (2.6) we find 

n 

[h* I kid. I ‘Max II A' (x) II *t] 

Ir— 1 ^ 

(2.8) llexp[^t]n < e , 0<t<T. 


The above estimate serves as a discrete analogue to the standard energy 
estimate one has in the differential case, whose abstract version amounts to 

II V 2 (a(x)D -D A(x) ) II *T V 2 'Maxll A' (x) H *T 

(2.9)llexp[A(x)D^t]n < e ■ < e , 0<t<T. 


In the case of a non-constant A, the two estimates, (2.8) and (2.9), differ, 

n 

however, in the term h* appearing in the first; to guarantee 

k=l 

stability we assume this term to be bounded 


n 

(L) h* 1 kjcl 1 < Const. 

k=l ^ 


Condition (I.) assures us that the differencing operator D is in a sense 
local , thus reflecting the local nature of differentiation D^. 


The antisymmetry, periodicity (= circulant form) and a locality 
characterization are essential properties shared by the differential operator, 
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which the discrete differencing operator D should retain as well, in order 
for an energy estimate (= stability) to be still valid under the discrete 
framework# 

Regarding the locality requirement, we consider for example the centered 
divided differences for fixed values of s; these operators are clearly 
local as they employ information extracted from a fixed number of neighboring 
gridvalues; this is also reflected in their matrix representation which 
has a finite width, w, defined as (see (1.5)) 


w(D) = Max {k|d ^O}. 
Kk<n 


We have w(jO=D 2 g) = s. Indeed, for such finite-width operators, 
Const .h*"^, and hence the locality condition (I.) is satisfied 


Idkl < 


” 2 
h* I k|d, I < Const. w (D), 

k=l 

yielding stable approximations. As s increases, however, becomes a 

full matrix which falls to satisfy the locality condition (L). That is not 
to say that the approximation becomes unstable, since the locality condition 
we have obtained is sufficient yet unnecessary for stability. Sufficient and 
necessary locality conditions which guarantee stability are, as much as we are 
aware, not known; we expect, however, that a locality condition requiring 
h*|d^| to decay faster than 1/k, k=l,***,n, is optimal for L 2 -stabillty of 
the general variable-coefficients problem. (We remark that in the constant 
coefficient case [and in general, with a definite coefficient A(x) , where 
multiplication first by A”^(x) will bring us back to the former case], we 


have MaxlIA'(x)ll = 0 and hence stability follows Independently of a locality 
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restriction, llexp[^t]l < 1, see (2.8).) To overcome the above difficulty, 
arising with "nonlocal" methods, two standard types of remedy can be employed 
— skew-symmetric differencing and introducing dissipation via smoothing; we 
discuss them next. 


3. Skew-Symmetric Differencing and Smoothing 

The spatial part of the differential system (0.1) is — apart from low 
order terms — a skew-selfajoint one 


A(x)D + B(x) = [a(x)D -H) A(x) ] + [b(x) -V2A'(x)]; 

X XX 


skew-symmetric differencing is based on exploiting this formalism. Rewriting 
(0.1) in the form 

3^u(x,t) “ { V2 [A(x)3^u(x,t)+9^(A(x)u(x,t) )] + [b(x)- (x) ] }u(x, t) ; 

and replacing the spatial derivative by its discrete counterpart, we end up 
with the approximation 

•= { V2 [aD4D^] + [b- V2i6']}v(t). 

The stability of the approximation in its skew-symmetric form is immediate: 
the first term inside the curly brackets is antisymmetric and hence its 
exponent is bounded by 1; by the perturbation lemma, therefore, the exponent 
of the sum of the two terms Inside the curly brackets is bounded by the 
exponent of the norm of the second 
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IIB- V2A'IIT 

llexp[{ V 2 [^+m] + [b- V 2 A' ]}t] II < e , 0<t<T; 


this is the exact energy estimate one has in the differential case. Thus 
skew-symmetric differencing, which is also available for a wide class of 
nonlinear problems, [17], maintains stability by retaining essential 
properties of the whole spatial operator A(x)D^ + B(x) rather than 
differentiation Itself; this is done, however, at the expense of doubling the 
total amount of work required. 

An alternative less expensive procedure to maintain stability is 
smoothing, a topic which the rest of this section is devoted to. We start by 
going back to estimate (2.6) where we were left with bounding the s)rmmetric 
part of Re(^) = V 2 (^-DA) • 

Employing the spectral representation of D, which we write as 
D = -gee (1.6), we obtain 

(3.1) I/2 (^-DA) = I/2 [a(n ^^"2 f)*A(n ^^^2 p) _ [n ^''2 f)*A(n ^^2 f)a] ; 


multiplying by N ^2 y qjj t-i^e left, by (n ^2p^ on the right, and observing 

1 /, 

that N is unitary (e.g. (A. 6) below), we find that the matrix above is 

unitarily similar and therefore equal in norm to 


(3.2) II V 2 (j^-DA) II = V 2 I|{[n^'^2p]a(n^^2p]*|a _ a{ (n ^^^2 p)a(n ^^^2 p)* } n . 


Next, we turn to examine the matrix inside the curly brackets, whose (p>q) 
block entry is given by 


p»q v=o 


(3.3) 
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using the Fourier expansion 


A(x) = I 


A(o)) 


1 f -io)5 

2¥ J ® 


OJs-oo 


A(?)dC, 


5=0 


it can also be expressed as 


(3.4) 


1 V 

n' ^ 
” v=0 


^i(q-p)vh 


i- I ( I A(«)e 


N 


v=0 a)=-<» 


a)=s-oo v=0 


I A(p-q+jN). 
j=-“ 


Having the representation of {(n ^2f)a(n | (3.4) and recalling the 

diagonal structure of A In (1.7), we conclude on account of (3.2) that the 
matrix V2 (AD-BA) is equal in norm to the matrix whose (p,q) block entry is 
given by 

^ ^ °° A 

(3.5a) ^-X^P ^ A(p-q+jN) 0<p,q<N-l. 

j=-“ 


We note that the locality condition (TL.) can be deduced again at this stage, 
if we are to proceed as follows: from (l*7b) we find 


(3.5b) 


/qO _ x(P') 


= 2i* ^ d • (sin(q^kh)-sin(p^kh) ) ( 
k=l ^ 


Since 

lsin(q"’kh)-sin(p''kh) 1 < k*h*Min[ 1 p-q I ,N- |p-q | ] , 


the matrix in (3.5) is bounded entrywise and therefore in norm, by the matrix 
whose (p>q) block entry is given by 
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2h 


I k|d^| •{Mln[lp-q| ,N-|p-qI ]• I HA(p-q+jN) H } 
k=l 


The matrix in the above curly brackets is a circulant one. As before, its 

norm does not exceed the absolute value sura of its elements along the first 

row (p-0) — see Corollary (A. 8) below; this sura in turn can be estimated in 

terras of the derivatives norm of A(x). Thus, assuming the locality condition 
n 

— h* 1 k|cL I < Const. — we conclude that V 2 and hence exp[^t], 

k=l 

0<t<T, have bounded norms, i.e., stability. 

The merit of the representation (3.5) lies, however, in the possibility 
of expressing a locality condition in terras of the amplification blocks 
associated with D, ^*1 , rather than its entries d, •! . To this end we 

proceed as follows: 

The matrix in (3.5) is written as the sum of two — the first takes the 
zero j-index which we rewrite as 


(3.6a) 


-(• 


q " p 


) 


•(p-q)*A(p-q); 


the second takes the rest of the j-lndices 

(3.6b) A(p-q+jN). 

j^O 

It is the property of the finite difference methods that the first matrix in 
(3.6a) is bounded. For, \ = 2i •ld^sin( j^kh) represents the discrete 

dif f erentation of the mode and as such, the order of magnitude of the 

difference should not exceed Const. | q'^-p^ I . Hence the matrix 

in (3.6 a) is bounded entry wise and therefore in norm, by the matrix whose 
(p>q) element is given by 
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Const. I p-ql • IIA(p-q) 0 ; 

the norm of such a Toeplltz matrix — see Corollary (A. 11) — does not exceed 
N-1 

Const.* I |u|llA(w)ll, which in turn, can be hounded by the norm of the 

( j )=0 

derivatives of A(x) . Regarding the boundedness of the second matrix in 

(3.6b), we note that for p-q bounded away from jN,j * 0, say |p-q| 

< 0N,0 < 1, we have H I A(p-q+jN) H < C and hence for these nonextreme 

j^O 

indices , the entries in (3 • 6b) are a' priori bounded — in fact , they are 
negligibly small. For the rest of the Indices, when Ip-ql N, i.e., when 
i p*' + q^ n, we must require the boundedness of \ Thus the 

locality condition amounts to the boundedness of the amplification 

blocks . ^ ^830C^-^ted with the high frequencies Ij^l n. If this is 

the case, the matrix V2 ^ts unitarily similar representation (3.5) 

is bounded and stability follows from (2.6) . 

The above situation is typical for all discrete methods, whose accuracy 
is determined by the exactness of differentiating the low modes, X'^^ ij^> 
while for their stability we need the boundedness of IX'^^ '| associated with 
the highest modes, |j^| The combination of the two guarantee 

convergence, as the low modes carrying most of the information are accurately 
represented, while the highest modes are inaccurately represented, yet 
stability assures us that they are not amplified and hence rapidly tend to 
zero, just as is the case in the differential problem* 

The two requirements — accuracy and stability — are well accommodated 
in difference methods having finite degree of accuracy; consider for example 


^^^It should be emphasized that this stability restriction is, of course, only 
sufficient* Its necessity is still an open question* 
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“I 

the second accurate method where we have = Ih sin(j^h), see (l«3b2)> 

and hence ^ ij' for |j^I *^0, i.e., accuracy, yet 

I A^J ^1 = |h"^sin(j^h) I < Const. for Ij^| ^ n, i.e., stability. The 

situation is less favorable, however, for highly accurate raehtods (of order 

N or more); the accuracy requirement A^*^ ^ ^ ij^ for the highest modes 

(i 

contradicts the stability restriction < Const. as originated from 

the locality conditon. Observe that this latter contradiction still leads to 
a bound proportional at most to N, which corresponds in the differential case 
to the familiar situation of "losing one derivatlve."^^^ 

The smoothing procedure aiming at dissolving this contradiction by 
bounding the amplification factors associated with the high frequencies (or 
more generally — the modes which these amplification factors multiply), yet 
leaving the lower accurate modes unharmed. For example, consider the so 
called Shuman filtering where 


is applied to the right hand side of (2.1a). In the Fourier space, it amounts 
to the further multlpllction of the i' mode, by V 2 * (l+cos( j"’h) ] ; that 

is 

Vj^ — ^^2 (l+cos(j ^h) )*Vj^. 


In other words, our smoothed discrete dif f erentatlon operator 
the form 


'^Shuman 


takes 


(3) 


In fact, as we shall see later on, we have a loss of "one-half" derivative 
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D_ = NF An_. F 

~Shuman '^Shuman 

with 

^Shuman “ (l+cos(-nh) ) • • • , V 2 (l+cos(((N-l-n)h) 

/ 4 -) 

which merely says that the amplification factors 1 were replaced by 

) . I /2 [i+cos(j'’h) ) . For the highest modes we now have the desired 
boundedness — in fact 1 ^ * V 2 (l+cos( j"h) ) | ~ 0 for |j"| ~ n. This is 
done, however, at the expense lowering the overall accuracy to a second one 
— ^ • V 2 (l+cos( j"h) ) “ Ij" +^[h^] for Ij"! ~ 0. In general, a 

linearly smoothed discrete dlfferentation operator may take the form 

(3.7a) D* = HF*A^ F, h* = Mi 

with 

(3.7b) dlag[a^““^ ‘I^, 

The requirement of both accuracy and stability can be now put in the concise 


form 




(3.8) ^ = } 

for 

ij'i 

1 bounded away from n (=* accuracy) 

(+0 

for 

ij'i 

1 + n (= stability) 

In [12], Majda et. al. 

advocated 

the use of exponential cut-off 


c>- smoothing when dealing with the propagation of singularities in linear 
problems# In [11], Kreiss and Oliger suggested a nonlinear smoothing, whose 
linearized version amounts to a polynomial cut-off of degree >2# In fact, a 
polynomial cut-off of degree one or more will suffice to compensate for the 
loss of one derivative we have observed earlier# To work out this last case 


in some detail, fix 9 < 1 and let 



(3.9) 


,cn 


Const.(| j^I-0n) 


-1 


in < On 

®n < Ij"! < n. 


(I"*) (1") d"') 

The adjusted amplification factors are now given by A'J ' — ► A'-’ 'o'"* A 
fixed portion of the N frequencies Is left unchanged maintaining the 

original order of accuracy. Regarding stability, we refer back to the real 

symmetric part ^ In Its unltarlly equivalent form (3.5) , which la written 
as the sum of two, see (3.6): the first 


- )-(p-q)-A(p-q) 


q -p 


is bounded by the norm of the derivatives of A(x) as we argued before; the 
second matrix 

3*0 

I Q 

is likewise bounded* Indeed, for Ip-q| < — 2 — entries are negligibly 
small — they are bounded by N* I lA(p-q+jN) * ^ C oN . For 

1+6 J*® 

Ip-ql > — - — *N we either have p > (l+9)n, q < (l-0)n, i*e. p^ > 0n, 

q^ < -0n, or, the roles of p and q are reversed# In either case 

Ip^l > 0n, Iq^I > 0n and therefore the latter matrix is essentially bounded 

entrywise and therefore in norm, by the matrix whose (p,q) entry is given by 

> ^ T q- ' i ^ en ~ Ip" I • 1 9^ l>0«- 

A direct calculation shows that this matrix is indeed bounded in terms of the 


norm of the derivatives of A(x) • 



Part II. The Fourier Method 


• The Fourier Differenclne Operator 


As before, we let v(x) be a 27r-perlodic m-dimensional vector-function, 
whose values v^ = v(x^) are assumed known at the gridpoints 
x^ =» Vh, h =» to simplify the notation we consider first the case of odd 
number of gridpoints, N = 2n+l, v =* 0,1, •••,2n. By Fourier differentiation we 
merely mean differentiation of the trigonometric interpolant of these 
gridvalues. That is, one construct the trigonometric interpolant 


(4.1a) 


v(x) = I V e 
u)=-n 


where the discrete Fourier coefficients v are given by, compare (1.8), 


(4.1b) 


^0, “ N* 

v=0 


-iuvh 




The Fourier differentiation then takes the form 


(4.2) 


-V. 

-^(Xy) = I loJVj^e 
u>=s-n 


The above procedure consists of the following three basic steps. First, 
transforming from the discrete space v = * * * * *^2n^ ^ into the Fourier 

A A A 

space of amplitudes v = (v • 


(4.3) 


V = Fv; 

A.# ' 


next, differentiation in the Fourier space takes place: 
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with A denoting the block diagonal matrix 

(4.4a) Ap = diag[-ln*I^,-i(n-l)*I^,**»,i(n-l)*I^,in*I^]; 

finally, the differentiated amplitudes A^y are transformed back into the 
discrete physical space: 

3plv] = F'HApv], f"^ = NF*. 


Added altogether, the Fourier differencing operator F amounts to 
multiplication by 

(4.4b) F - NF*A^F, 


which can be efficiently Implemented by two FFT's and N scalar 
multiplications requiring SNlogN operations • 

An explicit representation of the Fourier differencing matrix, F, can be 
obtained by differentiating the interpolant formula, cf., [22, Chapter X] 




2n 

•I v^K(x-x^), K(0 = 

v=0 

slnf (n+ V? )5 

2 sln(V 2 ?) * 

giving 




(4.5) 

1 

II 

1-) 

(-1)’^"^ 

0<j,k<2n. 

2sln[(k-j)x/(2n+l) J ni’ 


Thus, it falls into the category of antisymmetric block circulant matrices 
discussed above, see (1.5), 


(4.6) 


^ ,(F) (F) ^ (-1)^'*'^ 

'--Jjk [k-j] m I 2sintJiTr/(2n+l)] 
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wlth a spectral representation given by (4.4). Indeed, a straightforward 
calcualtlon, cf., Forenberg [3], shows 


~F 


llmA» ; 
^~2s’ 

S-Voo 


that is , the Fourier differencing can be viewed as a special centered finite 

differencing^ based on an ever increasing number of gridpoints extended in a 

periodic way , F » Continuing with this point of view, we conclude 

3 >00 

that while the Fourier differencing enjoys an "infinite order of accuracy" — 
a statement to be made precise below — it is a nonlocal one. We would like 
to examine the role these properties play in the Fourier method, based on 
replacing spatial derivatives by Fourier differencing. We start by discussing 
the all Important aliasing phenomenon. 


5. Aliasing 

Let w(x) be a smooth 27r-perlodic m-dimensional vector-function, with a 
formal Fourier expansion 

(5.1a) w(x) = I w(o))e^‘^ 

0>=_oo 

where the Fourier coefficients w(w) are given by 
(5.1b) w(aj) = / w(5)e 

Its interpolant w(x) based on the sampled gridvalues w(Xy),v = 0,l,***,2n, 
is given by 

w(x) =• I Wj^e^“^ 
o^-n 


(5.2a) 
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wlth the discrete Fourier coefficients 

2x1 

(5.2b) ''oj “ N* ^ w(x lw|<n. 

v=0 

The relation between the Fourier coefficients w(w) of w(x) and the 

A 

coefficients of its interpolant w(x) , is contained in the following 


Aliasing Lemma For w(x) as above we have 


(5.3) 


yv <» A 

w = I w(aH-kN). 
k=-“ 


Proof Inserting (5.1a) into (5.2b) we obtain 


^ v=0 v=0 


By the assumed smoothness of w(x), summation can be interchanged, yielding 


w =» w(p)*~* I w]h ^ ^ w(w+kN), 

lj=«oo v=0 


as the second sum in the middle term is nonvanishing only for those indices 

U such that [y-tj] = 0, i.e., y = w + kN. This completes the proof. 

Next, we consider the error between the grldf unction w(x) and its 

equidistant Interpolant w(x) . Rewriting w(x) = [!+)! ]w(03)e^^^, and, 

Io)|<n |w|>n 

with the help of the aliasing lemma. 


w(x) = I w(o))e^ ^ + I ( I w(aH-kN)]e^ 
|w|<n |w|<n k^O 
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we see that the difference w(x) - w(x) Is given as the sum of two basic 
contributions: the first, the truncation error , consisting of the higher 

truncated modes for Iw| > n 

(5.4a) Truncation [w(x)j * w(w)e^^, 

Iw|>n 

and the second, the aliasing error consisting of the higher aliased modes 
which were folded back on the lower ones, |w| < n, because of the finite 
resolution of our grid 

(5.4b) Aliasing [w(x)] = -I { I w[(i>fk(2n+l) ] 

|o)I<n k^O 

Observe that while the truncation error invloves modes higher than n, the 
aliasing error Involves modes less or equal to n; hence the two are 

orthogonal with respect to each other, and the size of the difference 
w(x) - w(x) is given by 

(5.5a) Ilw(x)-w(x) = IITruncatlon(w) + HAllaslng(w) H^. 

By Parseval's relation, the two squared terms on the right are given 
respectively by 

(5.5b) llTruncation(w) |w(w)|^ 

(w|>n 

(5.5c) HAllaslng(w) I I w[<i>fk(2n+l) ] | ^. 

Io)|<n k^O 

In both terms only the high amplitudes — those associated with modes higher 
than n - are being summed. Since these high amplitudes tend rapidly to 
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zero, i*e*, for smooth w(x) we have Iw(<a)) | < C^(l+|<^|) for any Y > 0, 

It follows that the two terras have the sarae error contribution of order 

C *N^ Y+1)^ Likewise we find that the derivative of w(x) , 3 w(x) , differs 
Y X 

from the differentiated interpolant, 3^w(x) by 

119 w(x)-9 w(x)B^ = I |(i)I?|w(u) 1^ + I |w|*l I w[oH-k(2n+l) ] 1^ 

!w|>n |w|<n k^O 

(-Y+2) 

which is of order C^N • As pointed out above, the Fourier differencing 

of w(x) is in fact the exact differentiation of the interpolant w(x)« We 
therefore conclude that the error we coramit by differentiating w(x) rather 
than w(x) is of the negligibly small order C^h for any 6 > 0. It is in 
this sense that we say the Fourier differencing has "infinite order accuracy.” 
Finally, we use the aliasing lemma to show the isometry between the 
discrete and continuous space functions. Precisely, consider the discrete 
space of grldfunctions (yQ> * * ~ ^^0* * * * * ^2n^^ equipped with the 

discrete inner product (*»*) 

(5.6a) (X,z) = h* I 

v=0 

as the discrete analogue of the space of 2 tt - per iodic vector functions, 

y(x),z(x) with 

27T ^ 

(5.6b) (y(x),z(x)) = / z (^)y(5)d5. 

0 

The above mentioned isometry now takes the concise form 


(5.7) 


(y(x),z(x)) = (x»z) 
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Indeed, consider the scalar Zir-periodlc function w(x) “ 2ir*z*(x)y(x) • While 

A 

the left-hand side is, by definition, w(oo=*0), the right hand side is, by 

definition, According to the aliasing lemma, the two differ by the sum 

of amplitudes associated with aliased modes higher than 2n, w[k(2n+l)]* 

k^O 

This sum is vanishing, however, since w(x) being a trigonometric polynomial 
of degree 2n at most, contains no modes higher than 2n* 


6* Stability of the Fourier Method 

In this section we study the stability of the Fourier method where 
spatial differentiation in (0.1) is carried out by Fourier differencing* 
According to the perturbation lemma we can safely neglect the low-order term 
assuming B » 0, and hence our approximation takes the form 

(6.1a) 

with the operator L given by 


(6. lb) L = A(x)D^. 

Indeed, the stability question as discussed above is relevant here, i*e«, the 

( ) 

unboundedness of the amplification blocks, see (4.4a), « ij^*I 

F m 

requires smoothing of the highest modes. In agreement with the nonlocallty of 

the method, see (4.6), h* \ k|d^^^I =^(l/h). The representation given below 

k-1 

la from a somewhat different point of view, and In fact. It Is the one that 
motivated our discussion In Section 3 above. 

Ikiltlplylng (6.1a) by hv* and summing over all grldpolnts we obtain 
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2n-l ^ 2n-l ^ ^ 

h* I V \v (t) = h- I V Lv = 

v=0 ^ v=0 ‘ V 


Taking real parts on both sides and making use of the isometry concluded in 
Section 5, we find 


( 6 . 2 ) 


dt 


llvli' 


2n-l 


= 2Re[h I v^3^v^(t)] = 2Re[(I^,v)] = 2Re[(Lv,v)] 


v=0 


The crucial step now, involves splitting the right hand side into the sum of 
two terms: the first consists of the exact differentiation 


(6.3a) 


2Re[(Lv,v)] = ([l+L*]v,v), 


the second consists of the deviation from the exact differentiation 

(6.3b) 2Re[(Lv-Lv,v) ] ; 

that is we have 

(6.4) 2Re[(Lv,v)] = 2Re[(Lv,v)] + 2Re[ [Lv-Lv,v] ] 

in complete analogy to the splitting of the matrix in (3.5) into (3.6a) and 
(3.6b) as we introduced before. 

^ 2 

That the first term in (6.3a) is bounded by Const .HvH is a property 
solely of the differential operator L, called semi-boundedness, which can be 
easily verified in our case by integration by parts. 


(6.5) 


I ( [l+L*]v,v) 1 < Const • • llvll^ , 


Const « • Max II A' (x) II 
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in complete agreement with the exponential behavior Indicated in (2.9). Thus 
we are left with estimating the second terra in (6.3b). It is exactly this 
term which measures by how much we deviate from the differential energy 
estimate whose abstract version quoted in (2.9). 

To this end we recall that the difference between w = Lv and its 
interpolant w = Lv consists of two basic contributions — the truncation 
error (5.4a) and the aliasing error (5.4b). The point to note here is that 
the truncation error being the sum of modes higher than n, is orthogonal to 
the n-degree interpolant v, and hence its contribution to the deviation term 
(6.3b) is completely suppressed. In other words, it is solely the aliasing 
error in the representation of the differential operator L — or what amounts 
to the same thing, of the coefficient matrix A(x) — which determines the 
stability of the discrete approximation (6.1). To see how it comes about one 
compute the amplitudes of Lv as the convolution sum 

n yv A 

(Lv] (tA)) = iq*A(w-q)v 
q=-n ^ 

hence the aliasing error is given by, see (5.4b) 

Allaslng[Lv] = { I I iq*A[o)-q+k(2n+l) ]v 

I w] <n I q| <n ^ 

Multiplying by v and making use of Parseval relation we find 

(Lv-Lv,v) = (Aliasing[Lv] ,v) = -i*I v *q* ^ A[p-q+k(2n+l) ] v ; 

|p|,|q|<nP k^O *1 


taking the S3rmmetrlc part we finally conclude that 
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(6.6) 2Re(Lv-Lv,v] = v*{(q-p)* ^ A[p-q+k(2n+l ) ] }v 

|p|,!q|<nP k^O ^ 

~ 2 

Our aim Is trying to estimate the right-hand side in terras of HvH — by so 
doing, then together with (6.5) we will end up with an energy estimate 

(6.7a) ^*'^**^ ^ Const •• Hvli 

whose integration assures us the continuous dependence of the solution on its 
initial data, i.e., stability, see (2.5) 

(6.7b) llv(t)ll^ = llv(t)ll^ < K(t)*llv(0)ll^ i K(t)*llv(0)ll^. 

^ 2 

To assert that the right-hand side of (6.6) does not exceed Const. Ilvl! = 
Const. I V I for all possible amplitudes v, . is, by definition, 

I I - CO CO 

col<n 

equivalent to assert the boundedness of the matrix whose (p>q) entry is 
given in the above curly brackets 

(6.8) = (q-p)’ I A[p-q+k(2n+l) ] , -n<p,q<n, 

k^O 

compare with (3.6b). The above terras represent the pure effect of aliasing on 
the coefficient matrix A(x) — in the constant coefficients case, for 

A 

example, no aliasing occur, A(to) =0, w ^ 0, so the terras in (6.8) and hence 
that in (6.3b) are vanishing which agrees with the earlier deduced stability 
in this constant coefficients case. Returning to the general variable 
coefficients case we first note — regarding the (p>q) entry in (6.8) — 
that for Ip-ql bounded away from 2n, |p-ql ^ 9*2n, 9 < 1, these entries are 
negligibly small, since by the smoothness of A(x) we have 



-34- 


II I A[p-q+k(2n+l)]ll < C gN“^; 
k^O 

however, when Ip-ql approaches 2n, that is, when pin and qi-n or vice 

versa, A[p-q+k(2n+l ) ] contains the lower modes of A(x) whose amplitudes 
k^O 

are of size ^(1), and hence these entries are of size ^(N=2n+1 ) — the 

matrix whose (p>q) entry is given in (6.8) is , therefore, unbounded, no 

matter how smooth A(x) Consider for example the case where A(x) 

consists of only one mode — the only nonzero entries in (6.8) are then the 

(p>q) ” (±n,+n) ones, given respectively by +2nA(ci)=+l) , vdiich amount to the 

unboundedness of (Putting it in a different way, we see that in constrast 

to local finite-difference methods, compare (2*7), Re(AF) = V 2 (AF-FA) is 

unbounded, no matter how smooth A(x) is; indeed up to unitary similarity — 

the latter differ from by the bounded terra in (6*3a))« 

Nevertheless, the above unboundedness does not necessarily imply 

instability, as much as it indicates the shortcomings of the above method of 

proving it. We observe that the difficulty arises when trying to estimate 

the (p>q) entries with pin and qi-n or vice versa, in either case, when 

!pl>lql n. These aliased entries interact with amplitudes associated with 

high modes (6*6), which are expected — the method is stable — 

to be of a negligible small size. That is, despite the unboundedness of 

in (6.8), we can still have the boundedness of the aliased term in (6.6) 

provided a^ priori Information on the decay rate of n is at our 

2 

disposal; it is well known, however, that the L 2 “norm II vH is too weak to 
derive such an a' priori information. 

With this an mind, smoothing may be viewed as a procedure aiming at 
giving us such a' priori Information about the size of the highest 
amplitudes. Consider for example the case where A(x) consists of fixed 
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number, say r modes; then smoothing by cutting a fixed number of modes — 
the last r ones, = 0, Io)|>n-r — will guarantee stability, as the term 

in (6.6) will vanish in this case. (In particualr, when r - 1, one only 

A 

needs to estimate the last amplitude v^. In the case of even number of 
gridpoints, such an estimate exists, since F being an even order 

antisymmetric matrix has a double zero eigenvalue; this leads to the 
H^-stability in this case derived in [7].^^^) In general, a milder smoothing 
as the linear cut-off introduced in (3.9) will suffice for stability. 

In closing this section, we remark that by rewriting (6.6) in the form 

2Re[lv-Lv,v) = i*I /l+|p | *v *{- • Y A[p-q+k(2n+l ) ] } •/I+TqT*v , 

|p|,|q|<n ^ /l+lq|/l+|pl k’^O ^ 

where the matrix in the last curly brackets is bounded, then together with 
(6.5) we are lead to the estimate 

(6.9) < Const. II Vll ^ 1 , , "vll^i, = I (l+|o)l |2^ 

H H '2 |o)|<n “ 

That is, there is a loss of "one-half" derivative. If some dissipation is 

2 

present in the system, e.g., with L = A(x)D^ + the gain of one derivative 
from the second spatial differentiation dominates, and we end up with the 
desired stability, e.g. [11]* 


(4) 


See the appendix for details. 



- 36 - 


Part III* Oalerkln Methods 


7. The Galerkln Procedure 

In this section we start discussing the Galerkin procedure, whose basic 
idea is to reduce our infinite dimensional differential problem by projecting 
it on a finite-dimensional subspace. Let the latter be spanned by a system of 
linearly independent 2 tt- periodic functions (x) , -n < p < n. 

To project (0.1), we seek for approximation of the form 

n .s 

(7.1) v(x,t) = I v(q,t)4> (x) 

q=-n ^ 

satisfying 

(7.2p) p=-n,*“,n. 

Inserting (7.1) into (7.2), we obtain for the vector of generalized Fourier 

A A A 

coefficients, v(t) = (v(-n, t) , * • • ,v(n, t) the following system of ordinary 
differential equations 

AK /Si 

(7.3a) M3^v(t) = Gv(t); 


here, M and G are (2n+l)- dimensional block matrices whose (p,q) entry 
Is given respectively by 


(7.3b) 


Mnn = 

pq q p m 


pq q p 


The stability of the resulting system is a direct consequence of the 
semi-boundedness of the differential operator L, cf . (6.5) — Integration by 
parts yields 

Re(Lw,w) = V 2 ( [l+L ]w,w) < Const • !lw II 

A 

for any 27r-periodic vector-function w(x) . Indeed, multiplying v(p,t) by 
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(7.2p), adding and taking real parts we find 

(7. A) V 2 ‘^llv(t) = Re(-|^,v) = Re(Lv,v) < Const. II v(t)ll^; 

Integration of the last Inequality gives us the usual stability estimate — 
compare (2.9). 

Unless chosen with care, the basis functions to an Ill- 

conditioned mass matrix, M, whose Inversion required In (7.3a) can be still 
found numerically disastrous. The most extensively studied choices of basis 
functions which avoid such situations are essentially two. The first uses 
local base functions Inducing sparse, well-behaved mass matrices, leading to 
flnlte-dlfference/flnlte-element like methods; the second uses global, 
orthonormal base functions like 

-n<p<n, 

where the mass matrix reduces to the identity, M = !• We continue by 
discussing the latter case. 

The expansion we seek in (7.1), amounts now to the truncated Fourier 
expansion, whose Fourier coefficients y(t) are determined by the Galerkin 
procedure 

A /N 

(7.4a) = Gv(t) 

where G Is given by 

(7.4b) [g]_„ = / A(x)e“^^^"‘*^*dx = Iq'A(p-q), -n<p,q<n. 

pq Q 

(As before, see (6.1b), we have neglected the lower order terra assuming L = 
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A(x)D^ ) • We remark that Implementation of the Galerkln procedure can be 
carried out fast « l»e., using 0 (NlogN) operations, provided the exact 
Fourier coefficients A((») , |(o|<n, are given. For, the procedure consists of 

two basic steps: first, dlfferentatlon which Is translated here to 

multiplication by the diagonal matrix taking place, 

requiring N » 2n+l operations; and next, multiplication by A(x) reflected 
as a convolution sum In the Fourier space Is In order, which requires 
multiplication by the Toeplltz matrix A(p-q) — Indeed, multiplication by a 
general Toeplltz matrix can be carried out fast when first Imbedded Into a 
clrculant one, see the appendix for details. 

To obtain the Fourier coefficients 

(7.5) A(p-q) = T 

^ 0 

one may use different quadrature rules approximating the integral on the 
right-hand side* This in turn leads to a whole variety of discrete Galerkin 
methods which include the Fourier method as a special case* 


8* Discretization 

The Fourier-Galerkin procedure in a component-wise form reads 

^ n ^ ^ 

(8.1a) 3j.v(p,t) •» I A(p-q) *lq*v(q,t) 

q=-n 

where A(w) Is the Fourier coefficient 

A(m) = J e ^*^*A(x)dx. 
x“0 


(8.1b) 
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To approximate the Integral in the right-hand side, we use the trapezoidal 

27T 

rule, based on the N = 2n+l equidistant points x^ = vh,h ■ — , 

A , N-1 

(8.2) A((o) ~ I k(x )e ; 

” v=0 

since A(x) is assummed periodic, the trapezoidal rule serves our purpose as 
any other high-order quadrature rule — in fact, it is infinitely order 
accurate” in the precise sense discussed in Section 5 above, cf* [2, Section 
2.9]. 

Introducing the approximation (8.2) into (8.1a), we find that the terra 
A(p-q) is replaced by, see (3.3) 

(8.3) A(p-q) I A(x = [nfAF*] ; 

thus, the above discretization result in a system of ordinary differential 
equation for the vector of unknown amplitudes, still denoted here by y, 

(8.4) ^^y(t) = NFAF 4 v(t). 

t r 

This is exactly the Fourier method for the discrete Fourier amplitudes 

AAA . 

X(t) = FV(t) « (v^^(t) , • • • ,v^(t) — multiplication by on the left 

brings it back into its familiar form in the physical space, see (4.4b) 

(8.5) 3^v(t) = A(NF*ApF)v(t) = ^(t) . 

That is, the exact spatial differentiation is carried out on the interpolant 
V, compare (6.1a). 
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To sunmarlze, we have seen that equidistant approximation of (8.3) based 
on N grldpolnts reduces the Fourier-Galerkln procedure into the Fourier 
method; the difference between the two lies exactly in the aliasing term 

A 

^ A(p-q+jN) — according to (3.4), this is the exact difference between the 
right and left-hand sides of (8.3). Since the Fourier-Galerkin procedure was 
shown to be stable, we thus shed a different light on the conclusion that 
stability of the Fourier method is solely determined by aliasing errors. To 
suppress the latter, one may either smooth or, alternatively, discretize the 
integral on (8.1b) using more than N grldpolnts. We turn now to discuss the 
details of the latter case. 

Let M = (1+G)N be the number of grldpolnts x,, = vh,h ■ 
v=«0, 1, • • *M-1, and use to trapezoidal rule to approximate 

(8.6) A(p-q) ~ A(x 

” v=0 

when Inserted into (8.1a), the resulting system is given by 

(8.7a) ^ ^ A(x )e^^‘^"P^^'^]*lq*v (t). 

P q=-n v-0 ^ 

Here, we adopt the notation of the discrete Fourier coefficients for the 
computed amplitudes, v^(t) , in the spirit of earlier agruraents. Observe that 
the matrix whose (p,q) entry, -n<p,q<n, is given in the last curly brackets, 
is not a clrculant anymore as in the Fourier case where M = N, yet its 
multiplication as a Toeplltz one can be carried out fast. To verify 
stability, we rewrite the (p,q) entry in the last curly brackets with the 
help of (3.4) 
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es. 7b) ^ [ I A(p-q+jM)]»iq*v (t). 

P q=-n j=-“ '* 

As usual, we break the second summation into two parts 

00 

/V ^ 

I A(p-q+jM) - A(p-q) + I A(p-q+jM) ; 
j=-CO j#0 

the first corresponds to the semi-bounded differential operator and can be 
estimated as before, while the second represent the pure effect of aliasing 
which in this case is completely controlled since by the smoothness of A(x) 
we have 

I llA(p-q+jM)II < C (£N)'^, Y>0, -n<p,q<n. 

j^O ^ 

Indeed, a second look in (8* 7b) reveals that the approximation (8*7) can be 
viewed as the standard Fourier method based on M modes, the last (1 +g) 
of which were cut off (in the notation of (3*8), we have =» 0 for 

(1+G)”^M < !jl < M ) — such smoothing guarantees stability. 
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A. On Toeplltz and Clrculant Matrices 

In this section we record some well-known Information about Toeplltz and 
clrculant matrices which proves useful within the discussion above 

A block Toeplltz matrix ^ consists of m-dlmenslonal block entries, 
the (j »k) of which depends only on Its distance from the main 
diagonal, j » 

(A.l) 


Thus, a Toeplltz matrix Is completely determined by a (2N-1)- dimensional 
vector _t = entries, -(N-1)<£<N-1, are being 

m-dlraenslonal blocks. 

If further, the vector _t Is defined on Its negative Indices as the 
periodic extension of the positive ones, t_^ = l.e., only 

depends on (k-j)[modN], then the matrix Is a block clrculant 

°"e,^,[^]jk “ ‘^(k-j)[modN] 





- 43 - 


(A.2) ^ 


*^0 

*^2 • • • *^-2 

* T^-1 

Sj -1 


® N-2 

^-2 


• 

• 

• 

• 


^2 

^2 



Cl C2 

* * *^-2 %-l 

^0 


Thus, a circulant matrix is completely determined by a N- dimensional vector 
_c = (^0 * * * * ’^-1 entries, 0 < ^ < N-1, being m-dlmensional blocks* 

The essential Ingredient in studying circulant matrices, is that they 
admit the spectral representation 

(A. 3) ^(c) = 

Here, F denotes the block Fourier matrix, compare (1*4) 

(A.4a) 0<j,k<N-l 

with the conventional notation 


(A. 4b) 


i' = £-n, n = integral part of N/2, 


and is a block diagonal matrix given by 


[Ajjj - 

^ £=0 


(A.5) 
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Verlflcation of (A. 3) is a straightforward one - the entry of the 

right-hand side of (A. 3) amounts to 






p=0 


£=0 


N 


M«.l 

i. J •j.lp'CA-J-Wh, 

^ £=0 p=0 


the second summation on the right is vanishing unless ^+J-k = 0[modN], i.e., 
unless ^ = (k-j) [modN] where 


^y^ ip^(^+j-k)h ^ _ f(y?] 

i=0%=0^ i^=(k-j)[raodN] ~ ‘^(k-j)[modN] = Jjk* 


Consideration of the block identity matrix as a circulant one, with 

c = (l ,0 ,***0 ], gives us from (A. 3) and (A. 5) that 
— mm m 

(A.6) = (n^/2f)*(n^/2f). 


that is, the matrix N ^^2 p is a unitary one. Since the spectrum and the L 2 - 
norra of a matrix are invariant under such unitary transformations, it follows 
from (A. 3) that for general circulant matrices, ^ , both are identical with 
those of block diagonal A^. In particular we have 


Lemma (A. 7) For a block circulant matrix ^(c^) we have 


II ^(c) II 


N-l ij)l^ 

Max II I e *c-ll. 
0<j<N-l £=0 


(A.7) 



-45- 


Proof The norm of a block diagonal matrix is given by the largest norm 
of its diagonal entries* Cosmetic reindexing of these diagonal entries in 
(A. 5) gives us (A. 7). 

As an immediate corollary we have 

Corollary (A* 8) The norm of a scalar circulant matrix does not exceed 
the absolute value sum of its elements along its first row. 

Proof In fact, from (A. 7) we have the more general 

N-1 

(A.8) ll^(c)ll < I llcJI. 

Jl=0 ^ 

The corollary is just a restatement of that last Inequality for the scalar 

case, where c„ = c„ *I . 

Z Z m 

Next, we employ the information just obtained for circulant matrices, for 
Toeplitz ones, with the help of the basic 

Lemma (A. 9) Any N-dimensional block Toeplitz matrix can be imbedded into 
a 2N-diraensional block circulant one. 

Proof consider the block Toeplitz raattix 3^ = ^(_t) with 

^ and define the associated Toeplitz matrix 

where s can be any fixed block. It is readily 


verified that 
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(A.9a) 


is a 2N-diraenslonal block clrculant 








(A. 9b) 




in entrywlse form we have 


€ 


•-1 


'2-N 


N-1 


N-2 N-1 


1-N ^2-N 


'1-N 


N-1 


'1-N 


N-2 N-1 


'-2 -1 


■-2 


■-1 


■-2 -1 "0 


--2 -1 


'1-N 2-N 


• • • 


N-1 


N-2 


'1-N 2-N 


■•l-N 

s 

Si-1 

*Sl-2 


'-1 0 


Remark Rewriting ^ in (A. 9c) as ^”(_t ) , _t = ) clarifies that 

the Imbedding was made psooible by the process of periodic doubling . 

Making use of Lemma (A. 9), we have 


Corollary (A. 10) Multiplication of an N-dimenslonal block Toeplltz 
matrix can be Implemented 'fast', i.e., using ^(NlogN) block operations. 
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Proof We want to compute z =» where ^ Is an N-dlmensional 

Toeplltz matrix and w a given N-dimenslonal vector* For that purpose, 

and compute z* = w* = (w,0^)" - 

as being a clrculant, this last multiplication can be Implemented fast 

using its spectral representation (A. 3) with two FFT^s requiring ^(NlogN) 
operations* The first N components of z^ are then the desired vector z* 

Corollary (A* 11) For a block Toeplltz matrix ^(jt) we have 

N-l 

(A. 11) II II < Max II e 

0<j<2N-l Jl=-(N-1) 

Proof Imbed ) into ^(_c) with ^ ) we then have 

from Lemma (A* 7) 

2N-1 

ll^(_t)ll < ll^f(^)ll = Max II I e *c»ll; 

0<j<2N-l £=0 

Insertion of the specific values of the blocks c^ in this case, shows that 
the upper-bound on the right-hand side equals 

N-l 2N-1 ij£^ N-l lj£^ 

Max II I e *tj| + (-1) *0 + I e om™ “ ^ ® 

0<j<2N-l £=0 £=N+1 0<j<2N-l £=-(N-l) 

Remark Making use of the freedom in choosing the block s along the 
main diagonal of the associated Toeplltz ^ (which was taken to be zero 
above) , we sirailarily get 

N-l ij£^ 

ll^(^)ll < inf[ Max ll(-l)^*s + I e *t.ll]. 

® 0<j<2N-l £=-(N-l) 


imbed &■ into «’ - ( f 

\(f S' 
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Corresponding to Corollary (A. 8) we have 

Corollary (A* 12) The norm of a scalar Toeplitz matrix does not exceed 
the absolute value sura of its elements along its first and last rows. 


B. The Fourier Method - The Case of Even Number of Gridpoints 

The Fourier method is usually irapleraented with an even number of 
gridpoints, N = 2n; to be exact, with N being an integral power of 2, in 
which case the Cooley-Tukey variants of FFT are optimal. Here we record the 
slightly different formulas governing this case. 

2tT _ 7T 

Assume v^ are known gridvalues at x^ = vh, h = — = — , v=0, 1, * • •2n-l. 
Their Fourier differencing araouonts to differentiation of their trigonometric 
interpolant 

(B.la) v(x) = I 

w=-n 


Here, the double prime denotes, as usual, halving the first and last terms, 
and the discrete Fourier coefficients v^ are given by 


(B.lb) 


6J 


, 2n-l 
= i. V V e 
v=0 


-iojvh 


An explicit representation of the Fourier differencing matrix F 

transforming v = ’ * * *^2n-l ” ^^x^lx ’*”’®x'^lx ^ 

0 2n— 1 

be obtained by differentiating the Interpolant formula, e.g. [22, Chapter X] 


1 2n-l 

(B.2) v(x) = I VyK(x-Xy) K(5) = 


giving 


sln(n5) 
2tg( V2?) 
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(B.3) - -(-D*'"J*cot((k-j)ii/2n)*I^, 0<j,k<2n-l. 

As a block clrculant matrix. It admits the spectral representation 

(B.Aa) F = NP*A F 

with 

(B.4b) A = dlag[0«I ,-l(n-l)*I ,•••,0*1 ,*‘*,l(n-l)*I ]. 

r m in m m 

Observe that zero is a double eigenvalue in this case - this is necessairly so 
as F being an antisymmetric even dimensional matrix, having the other 
complex eigenvalues in pairs. The left eigenvectors corresponding to the 
double zero eigenvalue are 

(B.5a) = 0, = (l ,I ,•••,! )' 

^ m m m 

(B.5b) = 0, = (l ,-I ,••*,1 ,-T 

m m m m 

asserting the exactness of the differentiation (B.3) for v(x) *« Const, and 
v(x) *= cos(nx) , respectively (compare [7, Lemma l.l]). 

The Fourier method for (0.1) with L = A(x)D^, is of the form 

(B.6) ^ 

Stability analysis in this case is similar to that Introduced in Section 7 for 
the case of odd number of grldpolnts. That is, to estimate the real symmetric 
part o£ ^Lv,vJ, see (6.2), we use the aliasing formula which still reads, see 


(5.3) 
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w = w(iiH-kN), N=2n 
k=-» 

leading us to an examination of the aliasing term, see (6«6) 

(B.7) 2Ke(Lv,v) “ i*I v [(q-p)» I A(p-q+2kn) ] v . 

|p|,|q|<nP k^O ^ 

In this case, however, we have a priori Information about the last discrete 

A 

Fourier coefficient v^. To see how it comes about, multiply (B.6) by J[ on 
the left, and rename the new variable w » Fv for which we find 

\s(t) - F^(t); 

next, multiplication by on the left and using (B.5b) we conclude 

that (z^^^)*w(t) - which, by definition, coincides with 

^ 2 n*“ 1 Ak ^ 

w, (t) = y w^.cos(nx,,) - remains constant in time, w. (t) ** w. (t*0) 
v=0 ±n in 

Thus, returning to the aliasing terra in (B*7), it is enough to sura only the 
first (n-1) modes 


2Re(Lw-Lw,w) “ i*I w [(q-p) I A(p-q+2kn) ]w • 

Ip! , Iql <n-l ^ k^O ^ 

In particular, if A(x) contains only one mode, the vanishing right-hand side 
results in the desired energy estimate for JJ = which in turn amounts to 
an stability in this case — see [7]. 
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